home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / synchrotron.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  6.9 KB  |  275 lines
  1. /* specfunc/synchrotron.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_exp.h>
  26. #include <gsl/gsl_sf_pow_int.h>
  27. #include <gsl/gsl_sf_synchrotron.h>
  28.  
  29. #include "error.h"
  30.  
  31. #include "chebyshev.h"
  32. #include "cheb_eval.c"
  33.  
  34. static double synchrotron1_data[13] = {
  35.   30.364682982501076273,
  36.   17.079395277408394574,
  37.    4.560132133545072889,
  38.    0.549281246730419979,
  39.    0.372976075069301172e-01,
  40.    0.161362430201041242e-02,
  41.    0.481916772120371e-04,
  42.    0.10512425288938e-05,
  43.    0.174638504670e-07,
  44.    0.22815486544e-09,
  45.    0.240443082e-11,
  46.    0.2086588e-13,
  47.    0.15167e-15
  48. };
  49. static cheb_series synchrotron1_cs = {
  50.   synchrotron1_data,
  51.   12,
  52.   -1.0, 1.0,
  53.   9
  54. };
  55.  
  56. static double synchrotron2_data[12] = {
  57.   0.4490721623532660844,
  58.   0.898353677994187218e-01,
  59.   0.81044573772151290e-02,
  60.   0.4261716991089162e-03,
  61.   0.147609631270746e-04,
  62.   0.3628633615300e-06,
  63.   0.66634807498e-08,
  64.   0.949077166e-10,
  65.   0.1079125e-11,
  66.   0.10022e-13,
  67.   0.77e-16,
  68.   0.5e-18
  69. };
  70. static cheb_series synchrotron2_cs = {
  71.   synchrotron2_data,
  72.   11,
  73.   -1.0, 1.0,
  74.   7
  75. };
  76.  
  77. static double synchrotron1a_data[23] = {
  78.   2.1329305161355000985,
  79.   0.741352864954200240e-01,
  80.   0.86968099909964198e-02,
  81.   0.11703826248775692e-02,
  82.   0.1645105798619192e-03,
  83.   0.240201021420640e-04,
  84.   0.35827756389389e-05,
  85.   0.5447747626984e-06,
  86.   0.838802856196e-07,
  87.   0.13069882684e-07,
  88.   0.2053099071e-08,
  89.   0.325187537e-09,
  90.   0.517914041e-10,
  91.   0.83002988e-11,
  92.   0.13352728e-11,
  93.   0.2159150e-12,
  94.   0.349967e-13,
  95.   0.56994e-14,
  96.   0.9291e-15,
  97.   0.152e-15,
  98.   0.249e-16,
  99.   0.41e-17,
  100.   0.7e-18
  101. };
  102. static cheb_series synchrotron1a_cs = {
  103.   synchrotron1a_data,
  104.   22,
  105.   -1.0, 1.0,
  106.   11
  107. };
  108.  
  109. static double synchrotron21_data[13] = {
  110.   38.617839923843085480,
  111.   23.037715594963734597,
  112.   5.3802499868335705968,
  113.   0.6156793806995710776,
  114.   0.406688004668895584e-01,
  115.   0.17296274552648414e-02,
  116.   0.51061258836577e-04,
  117.   0.110459595022e-05,
  118.   0.18235530206e-07,
  119.   0.2370769803e-09,
  120.   0.24887296e-11,
  121.   0.21529e-13,
  122.   0.156e-15
  123. };
  124. static cheb_series synchrotron21_cs = {
  125.   synchrotron21_data,
  126.   12,
  127.   -1.0, 1.0,
  128.   9
  129. };
  130.  
  131. static double synchrotron22_data[13] = {
  132.    7.9063148270660804288,
  133.    3.1353463612853425684,
  134.    0.4854879477453714538,
  135.    0.394816675827237234e-01,
  136.    0.19661622334808802e-02,
  137.    0.659078932293042e-04,
  138.    0.15857561349856e-05,
  139.    0.286865301123e-07,
  140.    0.4041202360e-09,
  141.    0.45568444e-11,
  142.    0.420459e-13,
  143.    0.3232e-15,
  144.    0.21e-17
  145. };
  146. static cheb_series synchrotron22_cs = {
  147.   synchrotron22_data,
  148.   12,
  149.   -1.0, 1.0,
  150.   8
  151. };
  152.  
  153. static double synchrotron2a_data[17] = {
  154.   2.020337094170713600,
  155.   0.10956237121807404e-01,
  156.   0.8542384730114676e-03,
  157.   0.723430242132822e-04,
  158.   0.63124427962699e-05,
  159.   0.5648193141174e-06,
  160.   0.512832480138e-07,
  161.   0.47196532914e-08,
  162.   0.4380744214e-09,
  163.   0.410268149e-10,
  164.   0.38623072e-11,
  165.   0.3661323e-12,
  166.   0.348023e-13,
  167.   0.33301e-14,
  168.   0.319e-15,
  169.   0.307e-16,
  170.   0.3e-17
  171. };
  172. static cheb_series synchrotron2a_cs = {
  173.   synchrotron2a_data,
  174.   16,
  175.   -1.0, 1.0,
  176.   8
  177. };
  178.  
  179.  
  180. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  181.  
  182. int gsl_sf_synchrotron_1_e(const double x, gsl_sf_result * result)
  183. {
  184.   /* CHECK_POINTER(result) */
  185.  
  186.   if(x < 0.0) {
  187.     DOMAIN_ERROR(result);
  188.   }
  189.   else if(x < 2.0*M_SQRT2 * GSL_SQRT_DBL_EPSILON) {
  190.     result->val = 2.14952824153447863671 * pow(x, 1.0/3.0);
  191.     result->err = GSL_DBL_EPSILON * result->val;
  192.     return GSL_SUCCESS;
  193.   }
  194.   else if(x <= 4.0) {
  195.     const double c0   = M_PI/M_SQRT3;
  196.     const double px   = pow(x,1.0/3.0);
  197.     const double px11 = gsl_sf_pow_int(px,11);
  198.     const double t = x*x/8.0 - 1.0;
  199.     gsl_sf_result result_c1;
  200.     gsl_sf_result result_c2;
  201.     cheb_eval_e(&synchrotron1_cs, t, &result_c1);
  202.     cheb_eval_e(&synchrotron2_cs, t, &result_c2);
  203.     result->val  = px * result_c1.val - px11 * result_c2.val - c0 * x;
  204.     result->err  = px * result_c1.err + px11 * result_c2.err + c0 * x * GSL_DBL_EPSILON;
  205.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  206.     return GSL_SUCCESS;
  207.   }
  208.   else if(x < -8.0*GSL_LOG_DBL_MIN/7.0) {
  209.     const double c0 = 0.2257913526447274323630976; /* log(sqrt(pi/2)) */
  210.     const double t = (12.0 - x) / (x + 4.0);
  211.     gsl_sf_result result_c1;
  212.     cheb_eval_e(&synchrotron1a_cs, t, &result_c1);
  213.     result->val = sqrt(x) * result_c1.val * exp(c0 - x);
  214.     result->err = 2.0 * GSL_DBL_EPSILON * result->val * (fabs(c0-x)+1.0);
  215.     return GSL_SUCCESS;
  216.   }
  217.   else {
  218.     UNDERFLOW_ERROR(result);
  219.   }
  220. }
  221.  
  222.  
  223. int gsl_sf_synchrotron_2_e(const double x, gsl_sf_result * result)
  224. {
  225.   /* CHECK_POINTER(result) */
  226.  
  227.   if(x < 0.0) {
  228.     DOMAIN_ERROR(result);
  229.   }
  230.   else if(x < 2.0*M_SQRT2*GSL_SQRT_DBL_EPSILON) {
  231.     result->val = 1.07476412076723931836 * pow(x, 1.0/3.0);
  232.     result->err = 2.0 * GSL_DBL_EPSILON * result->val;
  233.     return GSL_SUCCESS;
  234.   }
  235.   else if(x <= 4.0) {
  236.     const double px  = pow(x, 1.0/3.0);
  237.     const double px5 = gsl_sf_pow_int(px,5);
  238.     const double t = x*x/8.0 - 1.0;
  239.     gsl_sf_result cheb1;
  240.     gsl_sf_result cheb2;
  241.     cheb_eval_e(&synchrotron21_cs, t, &cheb1);
  242.     cheb_eval_e(&synchrotron22_cs, t, &cheb2);
  243.     result->val  = px * cheb1.val - px5 * cheb2.val;
  244.     result->err  = px * cheb1.err + px5 * cheb2.err;
  245.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  246.     return GSL_SUCCESS;
  247.   }
  248.   else if(x < -8.0*GSL_LOG_DBL_MIN/7.0) {
  249.     const double c0 = 0.22579135264472743236;   /* log(sqrt(pi/2)) */
  250.     const double t  = (10.0 - x) / (x + 2.0);
  251.     gsl_sf_result cheb1;
  252.     cheb_eval_e(&synchrotron2a_cs, t, &cheb1);
  253.     result->val = sqrt(x) * exp(c0-x) * cheb1.val;
  254.     result->err = GSL_DBL_EPSILON * result->val * (fabs(c0-x)+1.0);
  255.     return GSL_SUCCESS;
  256.   }
  257.   else {
  258.     UNDERFLOW_ERROR(result);
  259.   }
  260. }
  261.  
  262. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  263.  
  264. #include "eval.h"
  265.  
  266. double gsl_sf_synchrotron_1(const double x)
  267. {
  268.   EVAL_RESULT(gsl_sf_synchrotron_1_e(x, &result));
  269. }
  270.  
  271. double gsl_sf_synchrotron_2(const double x)
  272. {
  273.   EVAL_RESULT(gsl_sf_synchrotron_2_e(x, &result));
  274. }
  275.